function savings_EWM_cubic = generate_results_cubic_rule(bsimax,sample,wave,tcost,covariate,winter_prevuse,SMC,fullsample,savedir,fbase)
% This function calculates cost savings or energy reduction from EWM cubic
% rules using the cplex package. 

% Inputs: 
% (1) bsimax: number of bootstrap runs
% (2) sample: cross-sectional dataset merged using elec_wave367_cross.csv
% and propensity scores and renamed as data_estimation_pooled_ps
% (3) wave: indicates whether the output is wave-specific (=3,6,or 7) or
% for the pooled sample (=0)
% (4) tcost: private marginal cost of implementing the program per household
%  >0 represents cost savings; =0 represents kWh reduction 
% (5) covariate: indicates covariate for analysis: income, size, vintage,
% minimum of baseline consumption, maximum of baseline consumption, or
% standard deviation of consumption
% (6) winter_prevuse: indicates whether baseline consumption is calculated 
% as the mean of consumption in winter months (Jan and Feb) or as the mean
% of specified pre-treatment periods.
% (7) SMC: =1 use social marginal cost; =0 use retail electricity price
% (8) fullsample: =1 if full sample, =0 if estimation sample
% (9) savedir: directory to save the savings table 
% (10) fbase: string used to indicate the propensity score and baseline
% months specification for output table 

% Output: 
% (1) savings_EWM_cubic: a table summarizing the savings point estimate,
% share of households to treat, total savings for the entire state, and 95%
% confidence intervals for the point estimate 
%% Call CPLEX
addpath 'C:\Progra~1\IBM\ILOG\CPLEX_Studio129\cplex\matlab\x64_win64'

% Setup CPLEX optimization options
opt = cplexoptimset('cplex');
opt.parallel = 1;
opt.threads = 4;
opt.simplex.tolerances.feasibility = 1e-08;
opt.mip.tolerances.integrality = 1e-08;
opt.mip.strategy.variableselect = 3;
opt.mip.strategy.nodeselect = 2;
opt.mip.strategy.lbheur = 1;
opt.mip.limits.cutpasses = -1;
opt.display = 'off';

%% Gather inputs
[X,Y,D,n,ps,Yscale,x1_wave,prevuse_wave,Xscale] = generate_inputs_cubic_rule(sample,wave,tcost,covariate,SMC);
opower_paper_wave=sample.opower_paper_wave;

%% Estimate g using IPW
g = Y.*(D./ps - (1-D)./(1-ps)); % elements of W(G) - W(G_0)

%% Solve for rule using cplex 
% first BS permutation- original sample
bsperm = [1:n]';
gn = g(bsperm,:);
Xn = X(bsperm,:);

% number of variables (including the constant)
k = size(X,2);
bsi=0;

tic1= tic;
% Solve for optimal treatment rule
    % generate input for cplex package
    [f,f_n,Aineq_d,Aineq_i,bineq,lb,ub,ctype,Ind] = generate_cplex_inputs(bsperm,g,X,bsi,k,gn,Xn);
    % Cost minimization with treatment decreasing in pre treatment usage
    [sol_pd, v_pd] = cplexmilp(f,Aineq_d,bineq,[],[],[],[],[],lb,ub,ctype,[],opt); 
    % with treatment increasing in pre treatment usage
    [sol_pi, v_pi] = cplexmilp(f,Aineq_i,bineq,[],[],[],[],[],lb,ub,ctype,[],opt);
    
    % get beta and v 
    if (v_pd < v_pi)
            beta = sol_pd(1:k,:);
            v    = v_pd;
        else
            beta = sol_pi(1:k,:);
            v    = v_pi;
    end
    % get solved assignment in_Ghat and cost_saving
    in_Ghat=(X*beta>0);
    this_cost_saving=nanmean(g.*in_Ghat)*Yscale;

fprintf('cplex solution takes %.2f sec\n',toc(tic1))
    
v_pd_itr=0;
v_pi_itr=0;
v_nd_itr=0;
v_ni_itr=0;

%% Two-sided 95% confidence interval for the welfare gain/ cost savings
% Bootstrap for EWM CI
if bsimax>0
% initiate outputs
in_Ghat_itr=nan(n,bsimax);
v_itr=zeros(bsimax,1);
v_pd_itr=zeros(bsimax,1);
v_pi_itr=zeros(bsimax,1);
v_nd_itr=zeros(bsimax,1);
v_ni_itr=zeros(bsimax,1);
scaled_ATT_bs_itr = zeros(bsimax,1);
ATE_bs_itr = zeros(bsimax,1);
W_bs_Ghat_itr = zeros(bsimax,1);

g_hat_bs=zeros(n,bsimax);
cost_savings=zeros(bsimax,1);

% specification suffix 
if (tcost)
    if (SMC)
        s_speci = 'smc';
    else
        s_speci = 'pmc';
    end
else
    s_speci='kwh';
end
% shared name for each bootstrap batch
s_now_batch=datestr(now(),'yyyymmdd_HHMMSS');
% Create a saving directory for the i_BS files
 fpath_BS = ['.\intermediate_data\EWM_results\BS\cubic\',...
     s_now_batch,'_',covariate,'_',s_speci,'\'];
if ~exist(fpath_BS,'dir')
   mkdir(fpath_BS); 
end

sprintf('bootstrap:%s%s%d',s_now_batch,covariate,tcost)

parfor bsi = 1:bsimax 
    tic1=tic;
    rng(bsi);
    dt_start = now();
    % DRAW BOOTSTRAP SAMPLE WITH INDEXES bsperm
    bsperm = randi(n,[n 1]); % DRAW BOOTSTRAP INDEXES FOR THE NEXT SAMPLE
    
    %% save g_hat
    g_hat = generate_g_IPW(Y,D,bsperm);
    
    %% optimal rule 
    [f,f_n,Aineq_d,Aineq_i,bineq,lb,ub,ctype,Ind_bs] = generate_cplex_inputs(bsperm,g,X,bsi,k,gn,Xn);
    % Welfare maximization (cost minimization) with treatment decreasing in pre treatment usage
    [sol_pd, v_pd] = cplexmilp(f,Aineq_d,bineq,[],[],[],[],[],lb,ub,ctype,[],opt); 
    % with treatment increasing in post treatment usage
    [sol_pi, v_pi] = cplexmilp(f,Aineq_i,bineq,[],[],[],[],[],lb,ub,ctype,[],opt);
    
    %%
    % "Welfare minimization" for two-sided bootstrap C.I.
    [sol_nd, v_nd] = cplexmilp(f_n,Aineq_d,bineq,[],[],[],[],[],lb,ub,ctype,[],opt);
    [sol_ni, v_ni] = cplexmilp(f_n,Aineq_i,bineq,[],[],[],[],[],lb,ub,ctype,[],opt);
    
    if (v_pd < v_pi)
            beta_bs = sol_pd(1:k,:);
            v_bs    = v_pd;
        else
            beta_bs = sol_pi(1:k,:);
            v_bs    = v_pi;
    end
    in_Ghat_bs=(X*beta_bs>0);
    this_cost_saving_bs=nanmean(g_hat.*in_Ghat_bs)*Yscale;
    % timing
    dt_cmpt = now();
    
    %% output
    v_itr(bsi,1)=v_bs;
    v_pd_itr(bsi,1)=v_pd;
    v_pi_itr(bsi,1)=v_pi;
    v_nd_itr(bsi,1)=v_nd;
    v_ni_itr(bsi,1)=v_ni;
    cost_savings(bsi,1)= this_cost_saving_bs;
    in_Ghat_itr (:,bsi) = in_Ghat_bs;
    
    %% obtain att, ate, and savings from applying the original G_hat on bootstrap sample
    [scaled_ATT_bs,ATE_bs,W_bs_Ghat] = ewm_rct_test_cubic_input(D,Y,ps,beta,X,in_Ghat,g,Yscale,bsperm,n);

    scaled_ATT_bs_itr(bsi,1) = scaled_ATT_bs;
    ATE_bs_itr(bsi,1) = ATE_bs;
    W_bs_Ghat_itr(bsi,1) = W_bs_Ghat;
    
    %% save each BS result
    dt_outp = now();
    dt=[dt_start,dt_cmpt,dt_outp];
    save_BS([fpath_BS,sprintf('%s_BS_%d.mat',s_now_batch,bsi)],...
        bsi,s_now_batch,v_bs,v_pd,v_pi,v_nd,v_ni,...
        in_Ghat_bs,g_hat,scaled_ATT_bs,ATE_bs,W_bs_Ghat,bsperm,dt);
    % progress markers
    fprintf('Boot #%d took %2.3f sec\n',bsi,toc(tic1))

end

end

%% two-sided 95% confidence interval for the welfare gain/ cost savings
vhats_p= -min(v_pd_itr,v_pi_itr);
vhats_n= -min(v_nd_itr,v_ni_itr);
vhats_abs = max([vhats_p vhats_n],[],2);
 
%% Export rule parameters for graphs

switch winter_prevuse 
    case 1
        s_winter = '_winter';
    case 0
        s_winter = '';
end
if (tcost)>0 % cost_savings
    if (SMC)
        s_cost = 'SMC';
    else
        s_cost = 'PMC';
    end
elseif (tcost)==0
    s_cost='kwh';
end

if fullsample==0
    filename_coefs=sprintf('coef_cubic_%s_%s%s_wave%1.0f.mat',covariate,s_cost,s_winter,wave);
else
    filename_coefs=sprintf('coef_cubic_fullsample_%s_%s%s_wave%1.0f.mat',covariate,s_cost,s_winter,wave);
end

filename = sprintf('%s_%s',fbase,filename_coefs);

save(fullfile(savedir,filename),'beta',...
    'in_Ghat','x1_wave','prevuse_wave','g','covariate','tcost','Yscale','n','X','Xscale',...
    'v','vhats_n','vhats_p','scaled_ATT_bs_itr','ATE_bs_itr','W_bs_Ghat_itr','Ind');

%% Output for table 
percent=mean(in_Ghat);
savings=nanmean(g.*in_Ghat)*Yscale;
if (tcost)
    total_savings=savings*n*12;
else
    total_savings=savings*n*12/1000;
end
ci_lb=(v-prctile(vhats_abs,95))*Yscale/n;
ci_ub=(v+prctile(vhats_abs,95))*Yscale/n;

savings_EWM_cubic=nan(1,4);
savings_EWM_cubic=array2table(savings_EWM_cubic,'VariableNames',{'Covariates','Share\ treated','Net\ cost\ changes','Total\ cost\ changes'}); 
format short g
savings_EWM_cubic.(1) = {covariate};
savings_EWM_cubic.(2) =round(percent*100);
savings_EWM_cubic.(3) =round(savings,2);
savings_EWM_cubic.(4) =roundn(total_savings,3);

ci_lb=round(ci_lb,2);
ci_ub=round(ci_ub,2);

ci_table=nan(1,4);
ci_table=array2table(ci_table,'VariableNames',{'Covariates','Share\ treated','Net\ cost\ changes','Total\ cost\ changes'});
cilb_string=string(ci_lb);
ciub_string=string(ci_ub);
ci_for_table=strcat('(',cilb_string,',',ciub_string,')');
ci_table.(1)="";
ci_table.(2)="";
ci_table.(3)=ci_for_table;
ci_table.(4)="";
savings_EWM_cubic=[savings_EWM_cubic;ci_table];


end